---
title: "Replication Material for: How technological change affects regional voting patterns. Part 2: Main Analysis"
output: html_document
author: Nikolas Schöll and Thomas Kurer
editor_options: 
  chunk_output_type: inline
---

This script analyzes effect of robots on vote shares of different parties. We consider counties (Landkreis, n=324) and federal, state and EU election. If multiple election in the same year, we only consider one, preferring national over state over EU (order of  turnout)  which gives a total of 5115 region-election pairs.
We use a two way fixed effect model where we have a fixed effect for each region and a fixed effect for each election.
Standard errors are clustered on the county level.


# Load data
```{r 2_1, echo=FALSE}
options(stringsAsFactors=FALSE,
        scipen=999)
pacman::p_load(tidyverse, 
               magrittr,
               broom,
               knitr,
               tinytex,
               ragg,
               sf, # maps
               ggrepel, # value labels in ggplot
               lfe, # FE models
               stargazer) # output tables 

opts_chunk$set(echo=FALSE,
               message = F,
               warning = F)

df = read_rds("../data/sample.RDS")%>% 
  filter(east_west == "west") # ONLY WEST GERMANY

parties = c("greens", "linke",  "spd", "fdp",  "cdu_csu", "authoritarian_right")
party_labels = c("Grünen", "Linke",     "SPD", "FDP", "CDU/CSU" , "Authoritarian right")
party_labels_graphical = c("Grünen",  "Linke",   "SPD", "FDP", "CDU/CSU" , "Authoritarian\nRight")

# create results folders if missing
dir.create("../results/figures", recursive = TRUE)
dir.create("../results/tables", recursive = TRUE)
```

# Descriptive numbers
How many different elections covered? 50 in total
```{r 2_2}
df %>% 
  filter(year <= 2017,
         !is.na(election_type)) %>% 
  ungroup() %>% 
  distinct(year, election_type, election) %>% 
  count(election_type)
```

How many counties (Landkreise und kreisfreie Städte)?
```{r 2_3}
df$region %>% unique() %>% length()
```
How many election-year pairs (total n for political results)?
```{r 2_4}
df$spd %>% unique() %>% length()
```
How many county-year pairs (total n for labor market results)?
```{r 2_5}
df%>% 
  filter(!is.na(log_robots)) %>% 
  nrow()
```
# Data Section
## Figure 1: Evolution of manufacturing share, robot penetration and ICT
### 1A: Manufacturing share of employment

```{r 2_6}
# calculate manufacturing share separately for West and East Germany
manufacturing_shares_germany = 
  read_rds("../data/sample.RDS")  %>% 
  filter(year >= 2000) %>% 
  mutate(Country = case_when(east_west == "west" ~ "West Germany", 
                             east_west == "east" ~ "East Germany")) %>% 
  group_by(Country, year) %>% 
  summarise(manufacturing_share = sum(emp.siab_manufacturing) / sum(emp.siab) * 100)

# load manufacturing shares from other countries
manufacturing_shares_other_countries = 
  read_csv("../data/ILO/manufacturing share across countries.csv")  %>% 
  select(Country = ref_area.label, year = time, manufacturing_share = obs_value) %>% 
  filter(Country != "Germany")

manufacturing_shares =
  manufacturing_shares_other_countries %>% 
  rbind(manufacturing_shares_germany)  %>% 
  mutate(Country = factor(Country, levels = c("West Germany", "East Germany", "Spain", "France", "United States", "United Kingdom")))

manufacturing_shares %>% 
  ggplot(aes(x = year, y = manufacturing_share, color = Country, linetype = Country)) +
  geom_line(size = 1.2) +
  scale_y_continuous(labels = scales::unit_format(unit = "%", accuracy = 1)) +
  labs(x = "Year", y= "Manufactuing share of employment") + 
  theme(legend.position = "none")
  
ggsave(
  filename = "../results/figures/Figure_1a.pdf",
  width = 3.5,
  height = 3.5)
```

### 1B: Robots

```{r 2_7}
# calculate robots penetration separately for West and East Germany:
robots_germany = 
  read_rds("../data/sample.RDS") %>% 
  group_by(region) %>% 
  arrange(region, year) %>% 
  mutate(Country = case_when(east_west == "west" ~ "West Germany", 
                             east_west == "east" ~ "East Germany")) %>% 
  group_by(Country, year) %>% 
  summarise(robots = sum(robot_count) / sum(emp.siab) * 1000)

# load robot counts in other countries:

# This requires access to IFR data base
# robots_per_country = list()
# for (country_code in c("US", "UK", "ES", "FR")){
#   robots_per_country[[country_code]] =
#     read_csv2("../data/IFR/raw_data/IFR raw data.csv") %>%
#     filter(country == country_code,
#            industry == "000")  %>%
#     select(country_code_2d = country,
#            year = Year,
#            robot_count = `Operational stock`)%>%
#     mutate_at("year", as.character)  %>%
#   mutate(robot_count = robot_count %>% readr::parse_number())
# }
# robots_per_country %<>% bind_rows()
# robots_per_country %>% write_rds("../data/IFR/robots_per_country.RDS")
robots_per_country = read_rds("../data/IFR/robots_per_country.RDS")


# load employment data:
total_employment = read_csv("../data/ILO/total employment across countries.csv") %>% 
  mutate(country_code_2d = case_when(ref_area == "DEU" ~ "DE",
                                     ref_area == "FRA" ~ "FR",
                                     ref_area == "ESP" ~ "ES",
                                     ref_area == "GBR" ~ "UK",
                                     ref_area == "USA" ~ "US")) %>% 
  select(country_code_2d,
         country = ref_area.label,
         year = time,
         total_employment = obs_value)
  
# combine data on each part of Germany and other countries, calculate robots relative to employment:
robots_combined = robots_per_country %>%
  mutate(year = year %>% as.numeric()) %>% 
  left_join(total_employment) %>% 
  mutate(robots = robot_count / total_employment) %>% 
  select(Country = country, year, robots) %>% 
  rbind(robots_germany) %>% 
  filter(year >= 1994) %>% 
    mutate(Country = factor(Country, levels = c("West Germany", "East Germany", "Spain", "France", "United States", "United Kingdom"))) 

# plot:
robots_combined %>% 
  ggplot(aes(x = year, y = robots, color = Country, linetype = Country)) +
  geom_line(size = 1.2) +
  labs(x = "Year", y = "Robots per 1000 workers") + 
  theme(legend.position = c(0.2, 0.75))

ggsave(
  filename = "../results/figures/Figure_1b.pdf",
  width = 3.5, 
  height = 3.5)
  
```

### 1C: ICT

```{r 2_8}
# calculate ICT penetration separately for West and East Germany:
ICT_germany = 
  read_rds("../data/sample.RDS") %>%
  mutate(Country = case_when(east_west == "west" ~ "West Germany",
                             east_west == "east" ~ "East Germany")) %>%
  group_by(Country, year) %>%
 summarise(ICT = sum(ICT*emp.siab) / sum(emp.siab))

# load ICT data from other countries
ICT_other_countries = 
  read_rds("../data/ICT_country_comparison.rds") %>%
  mutate(Country = case_when(country == "ES" ~ "Spain",
                             country == "FR" ~ "France",
                             country == "US" ~ "United States",
                             country == "UK" ~ "United Kingdom")) %>% 
  ungroup() %>% 
  select(Country, year, ICT)

ICT_comparison =
  ICT_other_countries %>% 
  bind_rows(ICT_germany) %>% 
  filter(year >= 1995) %>% 
    mutate(Country = factor(Country, levels = c("West Germany", "East Germany", "Spain", "France", "United States", "United Kingdom"))) 

ICT_comparison %>% 
  ggplot(aes(x = year, y = ICT, color = Country, linetype = Country)) +
  geom_line(size = 1.2) +
  labs(x = "Year", y = "ICT capital stock per worker in €1000") + 
  theme(legend.position = "none")

ggsave(
  filename = "../results/figures/Figure_1c.pdf",
  width = 3.5, 
  height = 3.5)
```

## Figure 2: Regional distribution of new technologies


Generate label for top 5 counties in terms of robot intensity and ICT penetration
```{r 2_9}
# extract 5 counties with highest robot intensity
regions_most_robots_2017 = df %>% filter(year == 2017) %>% 
  arrange(-robots) %>% 
  select(region, region_code, robots) %>% 
  top_n(5) %>% 
  pull(region_code)

# extract 5 counties with highest ICT penetration
regions_highest_ict_2017 = df %>% filter(year == 2017) %>% 
  arrange(-ICT) %>% 
  select(region, region_code, ICT) %>% 
  top_n(5) %>% 
  pull(region_code)

#clean up region names for nicer display
df %<>% 
  mutate(region = case_when(region == "Muenchen" ~ "Munich (district)",
                            region == "Muenchen, Landeshauptstadt" ~ "Munich (city)",
                            region == "Hamburg, Freie und Hansestadt" ~ "Hamburg", 
                            TRUE ~ region))

# add to main data frame
df %<>% 
  mutate(region_label = ifelse(region_code %in% regions_most_robots_2017, region, NA) %>%
           str_remove_all(", Stadt") %>% 
           str_remove_all(", Landeshauptstadt"),
         region_label_ICT = ifelse(region_code %in% regions_highest_ict_2017, region, NA) %>% 
           str_remove_all(", Stadt") %>% 
           str_remove_all(", Landeshauptstadt"))

#remove intermediary data frames
remove(regions_highest_ict_2017,regions_most_robots_2017)
```

Prepare map
```{r 2_10}
theme_map <- function(...) {
  theme_minimal() +
  theme(
    text = element_text(color = "#22211d"),
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.major = element_line(color = "#ebebe5", size = 0.2),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#f5f5f2", color = NA), 
    panel.background = element_rect(fill = "#f5f5f2", color = NA), 
    legend.background = element_rect(fill = "#f5f5f2", color = NA),
    panel.border = element_blank(),
    ...
  )
}

map <- sf::st_read("../data/maps/shape_files_map_kreise/vg2500/vg2500_krs.shp", stringsAsFactors = F) %>% 
  mutate(region_code = RS %>% as.character() %>% as.numeric()) %>% 
  left_join(df %>% 
              filter(year == 2017),
            by = "region_code") %>% 
  filter(region_code < 11000)

map <- cbind(map, st_coordinates(st_centroid(map)))
```

### 2A: Robots
```{r 2_11}
map %>% 
  ggplot() + 
  geom_sf(aes(fill=robots, label = region_label)) + 
  scale_fill_gradient(low = "white", high = "black", trans = "log",
                      labels = scales::number_format(accuracy = 0.1)) +
   geom_label_repel(aes(X, Y, label.size=NA, 
                        label = region_label),
                    fill = "white",
                    color = 'black',
                    size  = 3,
           point.padding = 0.5) +
  theme_map() +
  labs(fill = "Robots per\n1000 workers", x = NULL, y = NULL) 
ggsave(filename = "../results/figures/Figure_2a.pdf")
```

### 2B: ICT
```{r 2_13}
map %>% 
  ggplot() + 
  geom_sf(aes(fill=ICT * 1000, label = region_label_ICT)) + 
  scale_fill_gradient(low = "white", high = "black",
                      labels = scales::unit_format(unit = "€")) +
   geom_label_repel(aes(X, Y, label.size=NA, 
                        label = region_label_ICT),
                    fill = "white",
                    color = 'black',
                    size  = 3,
           point.padding = 0.5) +
  theme_map() +
  labs(fill = "ICT capital stock\nper worker", x = NULL, y = NULL) 

ggsave(filename = "../results/figures/Figure_2b.pdf")
```


# Main results
## Figure 3: Region-level exposure to technological change and Party Vote Shares

### 3A: Robots
```{r 2_14}
results = list()
for (outcome in parties){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ log_robots | region +  election| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
    mutate(model = outcome,
           Specification = "Robots + Region FE + Year FE")
  
  results[[paste(outcome, "b")]] = df %>%
    felm(formula = get(outcome) ~ log_robots + ICT | region +  election| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
    mutate(model = outcome,
           Specification = "Robots + ICT + Region FE + Year FE")
}


results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = parties %>% rev(),
                       labels = party_labels_graphical %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip() +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0))+ 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))

  ggsave("../results/figures/Figure_3a.pdf", height = 10, width = 10,units = "cm")
```

Calculate within-region standard deviation in robot intensity: +35% 
```{r 2_15}
df %>%  
  group_by(region) %>% 
  mutate(robots_within = log_robots - mean(log_robots)) %>% 
  pull(robots_within) %>% 
  sd() %>% 
  exp()
```
### 3B: ICT
```{r 2_16}
results = list()
for (outcome in parties){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ ICT | region +  election| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Region FE + Year FE")
  
    results[[paste(outcome, 2)]] = df %>% 
    #filter(election_type != "state") %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  election| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Robots + Region FE + Year FE")
}


results %>% 
  map_df(~bind_rows(.)) %>%   
  mutate(Specification = factor(Specification,
                                levels = 
                       c("ICT + Region FE + Year FE",
                         "ICT + Robots + Region FE + Year FE") %>% rev())) %>% 
  ggplot(aes(x= factor(model,
                       levels = parties %>% rev(),
                       labels = party_labels_graphical %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity", 
             position = position_dodge(width=0.2),
             aes(shape = Specification),
             size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=0, position = position_dodge(width=0.2)) +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()+
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0)) + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))

  ggsave("../results/figures/Figure_3b.pdf", height = 10, width = 10,units = "cm")
```

Within region standard deviation in ICT capital stock per worker in 1000 euros: 523,43 EUR.
```{r 2_17}
df %>%  
  group_by(region) %>% 
  mutate(ICT = ICT * 1000, # transform unit to 1 euro
         ICT_within = ICT - mean(ICT)) %>% 
  pull(ICT_within) %>% 
  sd()
```
## Figure 4: Region-level exposure to robots and employment eﬀects

### 4A: Robots
```{r 2_18}
outcome_list = c(
  "employment_manufacturing_population_ratio", 
  "employment_non_manufacturing_population_ratio",
  "employment_population_ratio")

outcome_labels = c("Manufacturing",
                   "Non-\nManufacturing",
                   "Total")

results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ log_robots  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + ICT + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate)) +
  geom_point(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=.2, position = position_dodge(width=0.9)) +
  geom_hline(yintercept = 0, color = "grey") +
  #ylab("Effect of robots on employment to population ratios (in %)") +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() 

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  #ylab("Marginal effect of regional robot intensity on regional employment") +
  ylab(NULL) + 
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0)) + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))
  
  ggsave("../results/figures/Figure_4a.pdf", height = 10, width = 10,units = "cm")
```

### 4B: ICT
```{r 2_19}
results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ ICT  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Robots + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  mutate(Specification = factor(Specification,
                                levels = 
                       c("ICT + Region FE + Year FE",
                         "ICT + Robots + Region FE + Year FE") %>% rev())) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  #ylab("Marginal effect of regional ICT capital stock\non regional employment") +
  ylab(NULL)+
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0))  + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))

  ggsave("../results/figures/Figure_4b.pdf", height = 10, width = 10,units = "cm")
```


##  Figure 5: Technological change and Regional Task Composition

### 5A: Robots
```{r 2_20}
outcome_list =  c("task_cognitive_non_routine", 
                  "task_cognitive_routine", 
                  "task_manual_routine",
                  "task_manual_non_routine")
outcome_labels = c("Cognitive\nnon-routine", 
                   "Cognitive\nroutine", 
                   "Manual\nroutine", 
                   "Manual\nnon-routine")

results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ log_robots  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + ICT + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate)) +
  geom_point(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=.2, position = position_dodge(width=0.9)) +
  geom_hline(yintercept = 0, color = "grey") +
  #ylab("Effect of robots on employment to population ratios (in %)") +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() 

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  #ylab("Marginal effect of regional robot intensity on regional employment") +
  ylab(NULL) + 
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0)) + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))
  
  ggsave("../results/figures/Figure_5a.pdf", height = 10, width = 10,units = "cm")
```

### 5B: ICT
```{r 2_21}
results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ ICT  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Robots + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  mutate(Specification = factor(Specification,
                                levels = 
                       c("ICT + Region FE + Year FE",
                         "ICT + Robots + Region FE + Year FE") %>% rev())) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  ylab(NULL)+
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0))  + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))

  ggsave("../results/figures/Figure_5b.pdf", height = 10, width = 10,units = "cm")
```


##  Figure 6: Technological change and Regional Skill Requirements

### 6A: Robots
```{r 2_22}
outcome_list =  c("emp_hq_share", "emp_mq_share", "emp_lq_share")
outcome_labels = c("High-\nskilled", "Mid-\nskilled", "Low-\nskilled")

results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ log_robots  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + ICT + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate)) +
  geom_point(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=.2, position = position_dodge(width=0.9)) +
  geom_hline(yintercept = 0, color = "grey") +
  #ylab("Effect of robots on employment to population ratios (in %)") +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() 

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  #ylab("Marginal effect of regional robot intensity on regional employment") +
  ylab(NULL) + 
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0)) + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))
  
  ggsave("../results/figures/Figure_6a.pdf", height = 10, width = 10,units = "cm")
```

### 6B: ICT
```{r 2_23}
results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ ICT  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Robots + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  mutate(Specification = factor(Specification,
                                levels = 
                       c("ICT + Region FE + Year FE",
                         "ICT + Robots + Region FE + Year FE") %>% rev())) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  #ylab("Marginal effect of regional ICT capital stock\non regional employment") +
  ylab(NULL)+
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0))  + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))

  ggsave("../results/figures/Figure_6b.pdf", height = 10, width = 10,units = "cm")
```




## Figure 7:  Cross-sectional correlations of regional employment shares and party vote shares in 2017 Federal Elections
```{r 2_24}
super_mechanism_var_list =
  list(list(mechanism_label = "skill",
            mechanism_var_list = c("emp_hq_share", "emp_mq_share", "emp_lq_share"),
            mechanism_var_labels = c("High-skill", "Mid-skill", "Low-skill")),
       list(mechanism_label = "employment",
            mechanism_var_list = c("employment_non_manufacturing_population_ratio", "employment_manufacturing_population_ratio"),
            mechanism_var_labels = c("Non-Manufacturing", "Manufacturing")),
       list(mechanism_label = "main_task",
         mechanism_var_list =  c("task_cognitive_non_routine", 
                  "task_cognitive_routine", 
                  "task_manual_routine",
                  "task_manual_non_routine"),
            mechanism_var_labels = c("Cognitive\nnon-routine", 
                   "Cognitive\nroutine", 
                   "Manual\nroutine", 
                   "Manual\nnon-routine")),
       list(mechanism_label = "education",
            mechanism_var_list =  c("emp_edu_high_share", "emp_edu_low_share", "emp_edu_no_share"),
            mechanism_var_labels = c("A-levels\n(Abitur)", "High\nschool", "No High\nschool"))
)

for (i in 1:length(super_mechanism_var_list)){
results = list()
for(mechanism_var in super_mechanism_var_list[[i]]$mechanism_var_list){

for (outcome in parties){
  results[[paste(mechanism_var,outcome)]] = df %>% 
    filter(year == 2017) %>% 
    felm(formula = get(outcome) ~ get(mechanism_var)  | 0 | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("get")) %>% 
  mutate(model = outcome,
         mechanism_vars = mechanism_var)
}
}
g = results %>% 
  map_df(~bind_rows(.)) %>% 
  mutate(mechanism_vars = factor(mechanism_vars,
                                 levels = super_mechanism_var_list[[i]]$mechanism_var_list,
                                 labels = super_mechanism_var_list[[i]]$mechanism_var_labels)) %>% 
  ggplot(aes(x= factor(model,
                       levels = parties %>% rev(),
                       labels = party_labels_graphical %>% rev()), 
             y = estimate)) +
  geom_point(stat = "identity", position = position_dodge(width=0.9), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=0, position = position_dodge(width=0.9)) +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip() +
  facet_grid(cols = vars(mechanism_vars))

if(super_mechanism_var_list[[i]]$mechanism_label == "education"){
  g = g + facet_grid(cols = vars(mechanism_vars),
             scales = "free")
}
g %>% print()
ggsave(paste0("../results/figures/Figure_7_", super_mechanism_var_list[[i]]$mechanism_label, ".pdf"), height = 10, width = 16,units = "cm")
}

```


# Appendix: Robustness checks


## Figure A.1: Political Parties in the Two-Dimensional Space
```{r 2_24_2}
df_CHES = read_rds("../data/sample_CHES.RDS")

df_CHES %>% 
  ggplot(aes(x=dim_lr, y=dim_galtan, colour=party, label=label)) + 
  geom_point() + 
  geom_text_repel() +
  geom_vline(aes(xintercept = dim_lr_mean), linetype="dotted") +
  geom_hline(aes(yintercept = dim_galtan_mean), linetype="dotted") +
  xlab("Economic Left-Right Dimension") + 
  ylab("Postmaterialist-Authoritarian Dimension") +
  scale_colour_manual(breaks = df_CHES$party, 
                      values = (as.character(df_CHES$partycolor))) + 
  theme(legend.position = "none")

ggsave("../results/figures/Figure_A1.pdf", height=8, width=8)
```

## Figure A.2: Technological change and Regional Education Levels
### A.2A Robots
```{r 2_25}
outcome_list =  c("emp_edu_high_share", "emp_edu_low_share", "emp_edu_no_share")
outcome_labels = c("A-levels\n(Abitur)", "High\nschool", "No High\nschool")

results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ log_robots  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
  mutate(model = outcome,
         Specification = "Robots + ICT + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate)) +
  geom_point(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=.2, position = position_dodge(width=0.9)) +
  geom_hline(yintercept = 0, color = "grey") +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() 

results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  ylab(NULL) + 
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0)) + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))
  
  ggsave("../results/figures/Figure_A2a.pdf", height = 10, width = 10,units = "cm")
```

### A.2B ICT
```{r 2_26}
results = list()
for (outcome in outcome_list){
  results[[outcome]] = df %>% 
    felm(formula = get(outcome) ~ ICT  | region + year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Region FE + Year FE")
  
  results[[paste(outcome, 2)]] = df %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Robots + Region FE + Year FE")
}

results %>% 
  map_df(~bind_rows(.)) %>% 
  mutate(Specification = factor(Specification,
                                levels = 
                       c("ICT + Region FE + Year FE",
                         "ICT + Robots + Region FE + Year FE") %>% rev())) %>% 
  ggplot(aes(x= factor(model,
                       levels = outcome_list %>% rev(),
                       labels = outcome_labels %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  ylab(NULL)+
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()  +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0))  + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))

  ggsave("../results/figures/Figure_A2b.pdf", height = 10, width = 10,units = "cm")
```

## Figure A.3: Technological Change and Labor Market Composition

```{r 2_26_2}
df_meachanisms = read_rds("../data/sample_mechanisms.RDS")

#additional variables
cols = c("population_16_64", "log_population_16_64",  "emp_rate_55_64", "training",  "reintegration", "highskill_share", "abitur_share", "average_age", "median_age", "birth_rate", "share19_24", "share25_44", "share45_64", "net_inmove_total", "net_inmove_18_24", "net_inmove_0_17_30_49")

# run analogous regression and plot
# robots

results = list()
for (outcome in cols){
  results[[outcome]] = 
    df_meachanisms %>% 
    felm(formula = get(outcome) ~ log_robots | region +  year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
    mutate(model = outcome,
           Specification = "Robots + Region FE + Year FE")
  
  results[[paste(outcome, "b")]] = 
    df_meachanisms %>%
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("robots")) %>% 
    mutate(model = outcome,
           Specification = "Robots + ICT + Region FE + Year FE")
}


results %>% 
  map_df(~bind_rows(.)) %>% 
  ggplot(aes(x= factor(model,
                       levels = cols %>% rev(),
                       labels = cols %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity",
             position = position_dodge(width=0.2),
             aes(shape = Specification), size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width= 0, position = position_dodge(width=0.2)) +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip() +
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0))+ 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))
ggsave("../results/figures/Figure_A3a.pdf", height = 12, width = 18, units = "cm")

# ICT
results = list()
for (outcome in cols){
  results[[outcome]] = df_meachanisms %>% 
    felm(formula = get(outcome) ~ ICT | region +  year | 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Region FE + Year FE")
  
    results[[paste(outcome, 2)]] = df_meachanisms %>% 
    #filter(election_type != "state") %>% 
    felm(formula = get(outcome) ~ log_robots + ICT | region +  year| 0 | region) %>% 
    tidy() %>% 
    filter(term %>% str_detect("ICT")) %>% 
  mutate(model = outcome,
         Specification = "ICT + Robots + Region FE + Year FE")
}


results %>% 
  map_df(~bind_rows(.)) %>%   
  mutate(Specification = factor(Specification,
                                levels = 
                       c("ICT + Region FE + Year FE",
                         "ICT + Robots + Region FE + Year FE") %>% rev())) %>% 
  ggplot(aes(x= factor(model,
                       levels = cols %>% rev(),
                       labels = cols %>% rev()), 
             y = estimate,
             color = Specification)) +
  geom_point(stat = "identity", 
             position = position_dodge(width=0.2),
             aes(shape = Specification),
             size = 2) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=0, position = position_dodge(width=0.2)) +
  ylab(NULL) +
  xlab(NULL) + 
  scale_x_discrete() +
  coord_flip()+
  theme(legend.position = "bottom",
        legend.direction='vertical',
        axis.title = element_text(hjust = 0)) + 
      guides(color = guide_legend(reverse = TRUE),
             shape = guide_legend(reverse = TRUE))
ggsave("../results/figures/Figure_A3b.pdf", height = 12, width = 18, units = "cm")
```

## Prepare Appendix tables
```{r 2_27}
source("code_for_combined_tables.R") # load functions to combine tables

outcome_list = c(
  #party vote shares
  "greens", "linke",  "spd", "fdp",  "cdu_csu", "authoritarian_right",
  
  # labor market outcomes
  "employment_population_ratio", "employment_manufacturing_population_ratio", "employment_non_manufacturing_population_ratio"
)
```
## Tables A.1-A.6, A.13-A.15: Effect of robot intensity on party vote shares and larbor market outcomes
```{r 2_28, results='hide'}
table_counter = c(1:6,13:15)
for (i in 1:length(outcome_list)){ # loop over all parties and labor market outcomes
  outcome = outcome_list[i]
  print(outcome)
  
  #adjust time FE naming:
  if (outcome %in% parties){
    time_FE = "election"
    time_FE_label = "Election FE"
  } else{
    time_FE = "year"
    time_FE_label = "Year FE"
  }

  ### RUN REGRESSIONS ##########################################################
  ### OLS ###
  results = list()
  results[[1]] = df %>% felm(formula = get(outcome) ~ log_robots  | region + get(time_FE)| 0 | region) 
  results[[2]] = df %>% felm(formula = get(outcome) ~ log_robots + net_exports  | region + get(time_FE) |  0  | region ) 
  results[[3]] = df %>% felm(formula = get(outcome) ~ log_robots + ICT | region + get(time_FE) |  0  | region) 
  results[[4]] = df %>% felm(formula = get(outcome) ~ log_robots + gdp_per_capita  | region + get(time_FE) |  0  | region ) 
  results[[5]] = df %>% felm(formula = get(outcome) ~ log_robots + net_exports + ICT + gdp_per_capita | region + get(time_FE) |  0  | region ) 
  
  ### IV ###
  results_iv = list()
  results_iv[[1]] = df %>% felm(formula = get(outcome) ~  1 | region + get(time_FE)| (log_robots ~ log_robots_eu) | region) 
  results_iv[[2]] = df %>% felm(formula = get(outcome) ~  net_exports  | region + get(time_FE) |  (log_robots ~ log_robots_eu) | region ) 
  results_iv[[3]] = df %>% felm(formula = get(outcome) ~ ICT | region + get(time_FE) |  (log_robots ~ log_robots_eu)  | region) 
  results_iv[[4]] = df %>% felm(formula = get(outcome) ~ gdp_per_capita  | region + get(time_FE) |  (log_robots ~ log_robots_eu) | region ) 
  results_iv[[5]] = df %>% felm(formula = get(outcome) ~ net_exports + ICT + gdp_per_capita | region + get(time_FE) |  (log_robots ~ log_robots_eu)  | region ) 
  
  ### NON-LOGGED ROBOTS ###
    results_non_logged = list()
  results_non_logged[[1]] = df %>% felm(formula = get(outcome) ~ robots  | region + get(time_FE)| 0 | region) 
  results_non_logged[[2]] = df %>% felm(formula = get(outcome) ~ robots + net_exports  | region + get(time_FE) |  0  | region ) 
  results_non_logged[[3]] = df %>% felm(formula = get(outcome) ~ robots + ICT | region + get(time_FE) |  0  | region) 
  results_non_logged[[4]] = df %>% felm(formula = get(outcome) ~ robots + gdp_per_capita  | region + get(time_FE) |  0  | region ) 
  results_non_logged[[5]] = df %>% felm(formula = get(outcome) ~ robots + net_exports + ICT + gdp_per_capita | region + get(time_FE) |  0  | region ) 
  
  ### NON-LOGGED ROBOTS EXCLUDING OUTLIERS ###
  
  # identify outlier regions
  top_robots_regions =
    df %>%
    filter(year == 2017) %>% 
    slice_max(order_by = robots, n = 10) %>% 
    pull(region)
  
  results_non_logged_excl_outliers = list()
  
  results_non_logged_excl_outliers[[1]] = 
    df %>% 
    filter(!(region %in% top_robots_regions)) %>%
    felm(formula = get(outcome) ~ robots  | region + get(time_FE)| 0 | region) 
  
  results_non_logged_excl_outliers[[2]] = 
    df %>% 
    filter(!(region %in% top_robots_regions)) %>% 
    felm(formula = get(outcome) ~ robots + net_exports  | region + get(time_FE) |  0  | region ) 
  results_non_logged_excl_outliers[[3]] = 
    df %>% 
    filter(!(region %in% top_robots_regions)) %>% 
    felm(formula = get(outcome) ~ robots + ICT | region + get(time_FE) |  0  | region) 
  
  results_non_logged_excl_outliers[[4]] = 
    df %>% 
    filter(!(region %in% top_robots_regions)) %>% 
    felm(formula = get(outcome) ~ robots + gdp_per_capita  | region + get(time_FE) |  0  | region ) 
  
  results_non_logged_excl_outliers[[5]] = 
    df %>% 
    filter(!(region %in% top_robots_regions)) %>% 
    felm(formula = get(outcome) ~ robots + net_exports + ICT + gdp_per_capita | region + get(time_FE) |  0  | region ) 
    
  
  ### CREATE TABLE SEGMENTS ####################################################
  
  ### OLS ###
  command <-  paste('stargazer(results, covariate.labels = c("Robots", "Net Exports","ICT", "GDP per capita"), header=FALSE, out.header=FALSE, keep.stat=c("n","adj.rsq"), add.lines = list(c("Region FE", rep("X", length(results_iv))), c("', time_FE_label, '", rep("X", length(results_iv)))))')
  
  # store different segments of table separately
  begintable = stargazer.begintable()
  collabel = stargazer.keepcollabels()
  colnums = stargazer.keepcolnums()
  row_ols = stargazer.keeprow() # OLS row
  stats = stargazer.keepstats(command) #the stats part for the OLS (number of Obs, R2)

  ### 2SLS ###
  command = 'stargazer(results_iv , keep=c("robots"), covariate.labels=c("Robots"), header=FALSE, out.header=FALSE, keep.stat=c("n","adj.rsq"))'
  
  # store different segments of table separately
  row_iv = stargazer.keeprow() # IV row
  first_stage_fstat = paste0(c("First-stage F-stat", results_iv %>% map(~condfstat(.)) %>% unlist() %>% round(digits = 2)) %>% paste(collapse = " & "),"\\\\", collapse = "") # IV first stage F-stat

  ### NON_LOGGED ###
  command<-'stargazer(results_non_logged , keep=c("robots"), covariate.labels=c("Robots"), header=FALSE, out.header=FALSE, keep.stat=c("n","adj.rsq"))'
  
  # store different segments of table separately:
  row_non_logged = stargazer.keeprow()
  
  ### NON_LOGGED  EXCL OUTLIERS ###
  command<-'stargazer(results_non_logged_excl_outliers , keep=c("robots"), covariate.labels=c("Robots"), header=FALSE, out.header=FALSE, keep.stat=c("n","adj.rsq"))'
  
  # store different segments of table separately:
  row_non_logged_excl_outliers = stargazer.keeprow()
  FE_structure = 
    paste0(c("Region FE", rep("X", length(results))) %>% paste(collapse = " & "),
           "\\\\", "\n",
           c(time_FE_label, rep("X", length(results))) %>% paste(collapse = " & "),
           "\\\\") #Lines to indicate FE structure

  
  ### MERGE TABLE ##############################################################
  cat(begintable,
      collabel,
      colnums,
      "\\hline \\\\\\emph{OLS} \\\\", 
      row_ols, 
      "\\hline \\\\\\emph{2SLS} \\\\",
      row_iv,
      first_stage_fstat, 
      "\\hline \\\\\\emph{Non-logged robots} \\\\",
      row_non_logged,
      "\\hline \\\\\\emph{Non-logged robots exclude outliers} \\\\",
      row_non_logged_excl_outliers,
      "\\hline",
      FE_structure, 
      stats,
      "\\hline\\hline",
      "\\end{tabular}\n",
      file= paste0("../results/tables/Table_A", table_counter[i],".tex"), 
      sep="\n")
}
```

## Tables A7 - A12, A16-A18: Effect of ICT on party vote shares and labor market outcomes
```{r 2_29}
table_counter = c(7:12,16:18)
for (i in 1:length(outcome_list)){ # loop over all parties and labor market outcomes
  outcome = outcome_list[i]
  print(outcome)
  
  #adjust time FE naming
  if (outcome %in% parties){
    time_FE = "election"
    time_FE_label = "Election FE"
  } else{
    time_FE = "year"
    time_FE_label = "Year FE"
  }

  
  ### RUN REGRESSIONS ##########################################################
  ### OLS ###
  results = list()
  results[[1]] = df %>% felm(formula = get(outcome) ~ ICT  | region + get(time_FE)| 0 | region) 
  results[[2]] = df %>% felm(formula = get(outcome) ~ ICT + net_exports  | region + get(time_FE) |  0  | region ) 
  results[[3]] = df %>% felm(formula = get(outcome) ~ ICT + log_robots | region + get(time_FE) |  0  | region) 
  results[[4]] = df %>% felm(formula = get(outcome) ~ ICT + gdp_per_capita  | region + get(time_FE) |  0  | region ) 
  results[[5]] = df %>% felm(formula = get(outcome) ~ ICT + net_exports + log_robots + gdp_per_capita | region + get(time_FE) |  0  | region ) 
  
  ### IV ###
  results_iv = list()
  results_iv[[1]] = df %>% felm(formula = get(outcome) ~  1 | region + get(time_FE)| (ICT ~ ICT_eu) | region) 
  results_iv[[2]] = df %>% felm(formula = get(outcome) ~  net_exports  | region + get(time_FE) |  (ICT ~ ICT_eu) | region ) 
  results_iv[[3]] = df %>% felm(formula = get(outcome) ~ log_robots | region + get(time_FE) |  (ICT ~ ICT_eu)  | region) 
  results_iv[[4]] = df %>% felm(formula = get(outcome) ~ gdp_per_capita  | region + get(time_FE) |  (ICT ~ ICT_eu) | region ) 
  results_iv[[5]] = df %>% felm(formula = get(outcome) ~ net_exports + log_robots + gdp_per_capita | region + get(time_FE) |  (ICT ~ ICT_eu)  | region ) 
  
  
  ### CREATE TABLE SEGMENTS ####################################################
  ### OLS ###
  command <-  paste('stargazer(results, covariate.labels = c("ICT", "Net Exports","Robots", "GDP per capita"), header=FALSE, out.header=FALSE, keep.stat=c("n","adj.rsq"), add.lines = list(c("Region FE", rep("X", length(results_iv))), c("', time_FE_label, '", rep("X", length(results_iv)))))')
  
  # store different segments of table separately
  begintable = stargazer.begintable()
  collabel = stargazer.keepcollabels()
  colnums = stargazer.keepcolnums()
  row_ols = stargazer.keeprow() #OLS row
  stats = stargazer.keepstats() #the stats part for the OLS (number of Obs, R2)
  
  ### IV ###
  command = 'stargazer(results_iv , keep=c("ICT"), covariate.labels=c("ICT"), header=FALSE, out.header=FALSE, keep.stat=c("n","adj.rsq"))'
  
  # store different segments of table separately:
  row_iv = stargazer.keeprow() #IV row
  first_stage_fstat = 
    paste0(c("First-stage F-stat", 
             results_iv %>% map(~condfstat(.)) %>% unlist() %>% 
               round(digits = 2)) %>% paste(collapse = " & "),
           "\\\\", 
           collapse = "") 
  FE_structure = 
    paste0(c("Region FE", rep("X", length(results))) %>% paste(collapse = " & "),
           "\\\\", "\n",
           c(time_FE_label, rep("X", length(results))) %>% paste(collapse = " & "),
           "\\\\") # lines to indicate FE structure
  
  
  ### MERGE TABLE ##############################################################
  cat(begintable,
      collabel,
      colnums,
      "\\hline \\\\\\emph{OLS} \\\\", 
      row_ols, 
      "\\hline \\\\\\emph{2SLS} \\\\",
      row_iv,
      first_stage_fstat, 
  
      "\\hline",
      FE_structure, 
      stats,
      "\\hline\\hline",
      "\\end{tabular}\n",
      file= paste0("../results/tables/Table_A", table_counter[i],".tex"), 
      sep="\n")
}
```

Clean up
```{r}
rm(list = ls())
print("finished!")
```



